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Combining molecular dynamics and Monte Carlo simulation we study defect structures around an 
elongated colloidal particle embedded in a nematic liquid crystal host. By studying nematic ordering 
near the particle and the disclination core region we are able to examine the defect core structure 
and the difference between two simulation techniques. In addition, we also study the torque on a 
particle tilted with respect to the director, and modification of this torque when the particle is close 
to the cell wall. 
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I. INTRODUCTION 



Colloidal dispersions of small particles in nematic liq- 
uid crystals are a novel type of soft matter. Topological 
defects [0, |[ and additional long-range forces between the 
colloidal particles ^, ||] are immediate consequences of 
the orientational ordering of the liquid crystal molecules. 
The nematic-induced interparticle interaction brings a 
new range of effects to the system: supermolecular struc- 
tures 1^, , cellular structures |||, |j , and even a soft solid 
pof can be observed. Colloidal dispersions in liquid crys- 
tals also have a wide variety of potential applications [ pd] | . 

The subject of this paper is the liquid crystal ordering 
and equilibrium orientation of an elongated solid particle 
inside a uniformly aligned nematic liquid crystal. On the 
list of problems one has to clarify are: (1) the nematic 
ordering around the particle (including possible topologi- 
cal defects); (2) the type and strength of the orientational 
coupling between the particle and its aligned molecular 
environment; (3) the effect of confinement on the orienta- 
tional ordering of the particle, i.e. the equilibrium orien- 
tation of the particle close to the bounding surface. The 
solution to these problems is essential to understand the 
behavior of magnetic or nonmagnetic particles of colloidal 
size inherently present in biological liquid-crystalline tis- 
sues like cellular membranes. These problems also arise 
in ferroliquid crystals - the suspensions of single-domain 
ferroparticles in liquid crystals . 

The answer to the first question is known for spherical 
colloidal particles with homeotropic anchoringof the di- 
rector at the particle surface 0, § |, 0, M M' 0- 
Isolated particles provide a spherical confining geometry 
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for the liquid crystal. Sufficiently strong homeotropic 
anchoring induces a hedgehog defect with topological 
charge -1-1. The total topological charge of the whole 
system is zero, and an additional defect must be cre- 
ated to compensate the radial hedgehog. Two types of 
defect are possible: a hyperbolic hedgehog with a topo- 
logical charge —1, called a dipolar or satellite defect; or 
a — ^ strength disclination ring that encircles the spher- 
ical particle, called a quadrupolar or Saturn-ring defect. 
The dipolar (satellite) defect is a point defect, while the 
quadrupolar (Saturn ring) is a line defect. Theoretical 
and numerical work based on elastic theory ||, |l^, as 
well as computer simulation ^ , suggests that the 
dipole configuration is stable for micron-sized droplets. It 
is the one usually realized experimentally. The Saturn- 
ring defect appears if the droplet size is reduced or an 
external field is applied 

The same topological arguments are applicable for a 
non-spherical colloidal particle with homeotropic anchor- 
ing of the director at the particle surface. For an elon- 
gated particle with length L and transverse size D << L, 
and both L and D much greater than the dimensions of 
the molecules of the liquid crystal, one can have a discli- 
nation line of strength — 1 , a pair of disclination lines of 
strength — i, as well as the 'escaped radial' structure, in 
which the director bends over to become perpendicular 
to the particle surface [ p2[ . 

However, from the energetic point of view, the situa- 
tion is different from the case of the spherical particle. 
For the elongated particle, both defects are disclination 
lines. The elastic energy per unit length associated with 
a disclination of strength m is TiKm? ln(i?/ro), where R 
is the size of the sample and tq is a lower cutoff radius 
(the core size) |2^. This means that the free energy of 
a pair of — ^ disclinations is always smaller than that of 
a single —1 disclination. Therefore, one can expect that 
the pair of — ^ disclinations will always be a stable con- 
figuration. In principle, the —1 defect can still form a 
metastable state. 
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The answer to the second question is not known even 
in the framework of phenomenological (continuum) the- 
ory IHJ. The resuhs of the theory only indicate that, 
depending on the type and strength of anchoring, the 
equihbrium position of the particle may be either paral- 
lel or perpendicular to the liquid crystal director. The 
parameter governing the situation is the ratio of the par- 
ticle radius to the extrapolation length of the nematic 
liquid crystal. 

Obtaining an analytical expression for the torque (or 
elastic free energy) at arbitrary tilt angle 9 seems to be 
hardly possible due to the loss of symmetry of the direc- 
tor distribution and the presence of defects. Qualitative 
analysis shows that the proper argument for the free en- 
ergy should be (no • n)^ = cos^ 9 since the problem is 
bilinear in both the unperturbed director orientation no 
and the unit vector along the symmetry axis of the rod, 
n. A simple form of the free energy angular dependence 
has been proposed in Ref. p2| 



Ti_ + - T±_) cos^ 9, 



(1) 



which gives a sin 20 dependence for the torque and pre- 
dicts that the director response has a maximum at 9 = 
7r/4 and is absent a,t 9 — 0, 7t/2. However, it is clear that 
eqn (|^) is oversimplified. The defect structure changes 
while the particle rotates. The nematic ordering evolves 
in a complicated way that can hardly be approximated 
with a sin 29 dependence of the torque. At the same time, 
this dependence is vital for a macroscopic description of 
the system which treats coupling of colloidal particles 
with the nematic host via an effective potential. 

Recently, the effect of confinement on the orientation of 
an anisotropic colloidal particle has been predicted p4| . 
It has been shown that there is an 'entropic' torque on 
a hard rod-like particle dissolved in a solution of hard 
spheres when the rod is positioned close to the hard wall. 
The torque appears because of the density modulation of 
spheres near the wall and depletion forces between the 
wall and the hard rod. This torque might play an impor- 
tant role in the 'key-lock' principle in biological systems 
and provide an understanding of how a non-spherical 
'key' macromolecule can adjust its position and orienta- 
tion near the 'lock' macromolecule. We, therefore, expect 
the torque on the particle close to the wall to be different 
from the bulk-induced torque because of the additional, 
entropic contribution. 

In this paper, we present the results of molecular dy- 
namics (MD) and Monte Carlo (MC) simulations of the 
topological defects that appear in the nematic mesophase 
around an elongated colloidal particle. Using the MD 
technique, we also study the force and the torque on the 
particle suspended in the bulk of the nematic mesophase 
and the modification of this torque when the particle is 
close to the cell substrate. 

The paper is organized as follows. In section || we 
present the computational details and molecular models 
we use to simulate the liquid crystal mesophase and the 
interaction of the molecules with the particle surface and 
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FIG. 1: Studied geometry: a spherocylinder of length L 
and diameter (Jr is immersed in a liquid crystal host, which is 
modeled either as a solution of Gay-Berne particles (molecular 
dynamics) or 'spins' fixed on a cubic lattice (Monte Carlo sim- 
ulations). The symmetry axis of the spherocylinder is tilted 
with respe ct to the z axis. To study the defect structure 
(sec. [II A) we use a rod of infinite length positioned along 



the y axis, normal to the director. 



the cell substrates. 



Section III contains the results of 
density, director. 



the MD and MC simulations: density, director, order 
parameter maps, and order tensor profiles of the defects. 
Here we also present the results for the torque on the 
particle in the cell bulk and near the wall. Concluding 
remarks and comparison of the techniques are given in 
section IV. 



II. 



MOLECULAR MODEL AND SIMULATION 
METHODS 



Figure |^ shows the geometry of a single particle in the 
cell. The rod shape is chosen to be a spherocylinder, 
i.e. a cylinder of length L and diameter ar which has 
spherical caps of diameter ar- The particle is tilted in 
the zy plane. The orientation of the particle is specified 
by the angle 9 between the z axis and the symmetry axis 
of the cylinder. 



A. Molecular dynamics 

Molecular dynamics simulations were carried out using 
the soft repulsive potential, describing (approximately) 
ellipsoidal molecules 
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(2) 



Here gtj = {r.ij — Oij + (Jo)/oo; Vij is the center-center 
separation, (Tq ^ size parameter, Eq an energy parameter 
(both taken to be unity) and the orientation-dependent 
diameter aij is defined by 



<tI ~ 2 



1 -I- x(Mi • Mj) 



1 - x(m» • Mj) 
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where x = ('^^ ~ + 1); being the elongation. 

In this work we used k — 3 throughout. The orienta- 
tion dependence is written in terms of the direction of 
the center-center vector fij — rij/rij and the unit vec- 
tors Ui, Uj which specify the molecular symmetry axes. 
The potential (§) may be thought of as a variant of the 
standard Gay-Berne potential |2^, 26 1 with exponents 
^ = 0,1/ = 0. 

The systems consisted of = 64, 000 particles. A re- 
duced temperature k^T/eo — 1 was used throughout (for 
this model, the phase behavior is not sensitively depen- 
dent on temperature, as there are no attractive forces). 
The system size was chosen so that the number density 
of the liquid crystal far from the rod was paQ ~ 0.34. For 
this system, in the reduced units defined by ctq, eo and m, 
a timestep St = 0.004 was found suitable. The molecular 
moment of inertia was fixed as / = 2.5m<TQ. Periodic 
boundary conditions as well as slab geometry with walls 
confining the system in the z direction were considered. 

The interaction of molecule i with the rod was given 
by a shifted Lennard- Jones repulsion potential having ex- 
actly the same form as eqn (0), but with pij replaced by 
Pi = {\ri - rs\ - ar/2 + ao/2)/cro. Here = ^iU, where 
7i = sign.{n ■ Ti + Zr cos9) min {L/2,\n ■ ri + Zr coa9\), 
n = (0, sin 0, cos 6') is a unit vector along the symme- 
try axis of the rod, Zr is the distance from the center 
of the rod to the center of the coordinates. In slab 
geometry the interaction with the walls was given by 
the same formula, eqn (|^), replacing particle j by the 
wah w, setting = {\ziw\ - o-j^,/2 -I- ao/2)/ao and 
fj^uj = '^^ + (1 ~ - ef^). Ku, represents an effec- 

tive particle elongation as seen by the wall. We used 
— I which gives strong homeotropic orientation of 
the molecules at the wall ]l9[ . 

The radius and length of the rod were steadily in- 
creased from zero to the desired value during 10^ steps. 
Then the system was equilibrated for 10^ steps. During 
equilibration we scaled the velocities of the molecules to 
achieve k^T /ea = 1. 

The production run for every tilt angle of the rod was 
10^ steps. The force F and the torque M on the rod 
were calculated using the repulsive force f i from the rod 
on the particle i 

N 



N 

M = -E7.[/. xn]. 



(3) 
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B. Monte Carlo simulation 

Monte Carlo simulations were based on the Lebwohl- 
Lasher (LL) lattice model |2^. Within this model uni- 
axial nematic molecules (or, alternatively, close-packed 
molecular clusters [^) are represented by 'spins' fixed 
on a cubic lattice of spacing a. The spins are allowed to 



rotate freely, which reproduces the orientational behavior 
of the liquid crystal sufficiently well. 

To define the topology of the colloidal particle, a 
'jagged' cylinder of diameter Ur was carved from the cu- 
bic lattice, with its long axis fixed along the y axis of 
the coordinate system. The orientations of spins rep- 
resenting the particle were kept fixed during the simu- 
lation and were chosen in agreement with the desired 
boundary conditions at the particle surface, in our case 
homeotropic. At outer boundaries of the simulation box 
periodic boundary conditions were assumed. The total 
interaction energy for our model system consisting of ne- 
matic spins was calculated as 



U : 



E 

{i<j) 



P2(cos/3j. 



(4) 



Ui ■ Ui 



Here Uj 



with ^2(2^) — \{ix^ — 1) and cos/3i 
denotes the unit vector giving the orientation of the spin 
located at the ith lattice site. The sum in eqn (^ is taken 
over nearest neighbors only. The constants represent 
the interaction strengths and are denoted by e and for 
nematic-nematic and nematic-solid particle interactions 
respectively. 

The simulation box size was set to 30 a x 30 a x 30 a, 
which for the chosen cylinder diameter (tr^ — 10 a) 
amounts to 24600 nematic spins and to 840 spins repre- 
senting the surface of the solid particle. Our simulations 
started from a configuration with a random orientation 
of nematic spins. The final results did not depend on the 
choice of the starting configuration. 

The standard Metropolis scheme |^ was then em- 
ployed to update nematic spin orientations po| , |3l|| , main- 
taining a rejection ratio close to 0.5. The system was 
equilibrated during « 6 x lO'' MC cycles. After this 
6.6 X 10^ successive spin configurations were accumulated 
and used as input for the calculation of order tensor. 

In the simulation, temperature was set to T* = 
ksT/e = 1, which ensures the existence of the nematic 
phase (note that for a bulk sample the LL model exhibits 
a nematic-isotropic transition at T* ~ 1.1232 [Q). The 
strengths of nematic-nematic and nematic-solid particle 
interactions were set equal, ep = e, which corresponds 
to the strong anchoring regime with the extrapolation 
length of the order of a few lattice spacings a P3|. 



C. Order tensor 

For both techniques, the local order tensor Q{r) was 
calculated 



" k=l ^ ^ 



(5) 



where there are n molecules present in each bin. Sap is 
the Kronecker delta, (• • ■ ) denotes an ensemble average, 
a,l3 = X, y, z. Note that in the MC case we have n = 1 
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FIG. 2: MC simulation results: cross section of the director 
field n{x,z). The shading represents the value of the order 
parameter, S. A pair of — | defects has formed on the diago- 
nal. In the defect core molecules are (on the average) aligned 
in the xz plane; ordering is uniaxial with S < and the cor- 
responding eigenvector, n, is directed out-of- plane (along the 
long axis of the particle). 

(bins correspond to lattice points) and that averaging is 
performed over MC cycles only. Diagonalizing the Qafi 
tensor, for each bin, gives three eigenvalues, Qi, Q2, and 
Qs, plus the three corresponding eigenvectors. The eigen- 
value with the largest absolute value defines the order 
parameter S for each bin. The biaxiality P is then calcu- 
lated as the absolute value of the difference between the 
remaining two eigenvalues of the order tensor Q. 

III. SIMULATION RESULTS AND DISCUSSION 

A. Defect structure 

We started our simulations with a rod of infinite 
length, L = 00, positioned along the y axis, normal to the 
director. In this case the director rotates in the xz plane 
and we effectively have a two-dimensional situation. 

1. MC results 

Fig. H shows the director field and the order parameter 
map, in the plane perpendicular to the long axis of the 
colloidal particle. 

As concluded from topological considerations, either a 
— I strength disclination line or a pair of — i lines can 
form in the neighborhood of the particle. The —I line, 
however, does not seem to be stable and splits into a pair 
of — I lines during the MC evolution, even if it is taken as 
initial configuration in the simulation run. This behavior 
agrees with simple estimates of defect line free energies. 
Moreover, a stable 'escaped radial' structure was also not 
observed in our simulations. In MC simulations the pair 
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FIG. 3: MC simulation results: order tensor components Qi 
(triangles), Q2 (circles), and Q3 (squares) plotted across the 
defect. In the upper panel we plot the order parameter 5* 
(open circles) and biaxiality P (squares). The left-right asym- 
metry with respect to the defect core is due to the presence 
of the colloidal particle. 



of defect lines always forms close to one of the simula- 
tion box diagonals although the cross section of the col- 
loidal particle is axially symmetric (ignoring its jagged 
shape); see the director field shown in Fig. 0. This sym- 
metry breaking may be attributed to two effects of dif- 
ferent origin. The first one (and, according to our tests, 
the more important one) is the repulsion between defects 
maximizing the defect-to-defect distance (recall the pe- 
riodic boundary conditions), while the second one is a 
finite-size effect originating from collective fluctuations, 
resulting in a tendency to align the nematic along the 
simulation box diagonal Q. We believe, however, that 
these phenomena, as well as the presence of the colloidal 
particle, do not considerably affect any of the qualitative 
features characterizing the disclination line inner struc- 
ture. Moreover, the presence of the colloidal particle is 
reflected only in an enhancement of the degree of nematic 
order in the immediate surroundings of the particle. The 
inner structure of a defect line is further characterized 
by variations in order tensor components, Qi, Q2, and 
Qs, obtained after diagonalization of the order tensor Q, 
eqn (^. This is the most convenient way of describing 
the structure of the defect, because of the possible biaxi- 
ality and negative values of the uniaxial order parameter 
in the core region. 

Fig. § shows the Qi, Q2, and Q3 profiles plotted along 
the z axis through the left of the two disclinations shown 
in Fig. H. In Fig. || the disclination is located a,t d — 20a. 
Note that the left-right assymetry of the profiles with re- 
spect to the defect position is caused solely by the pres- 
ence of the colloidal particle. As shown by Fig. |^, the Qi 
component changes from its positive bulk value (« 0.6), 
coinciding with the value of the order parameter S) to 
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some negative value (« —0.3) after passing through the 
dischnation. At the same time, the Q2 component in- 
creases from a negative value (« —0.3) to a large positive 
value (~ 0-6), which roughly equals twice the absolute 
value of the negative one. This behavior is attributed 
to the director rotation by approximately 7r/2 when we 
cross the defect along the z axis; see Fig. |[ On the other 
hand, the value of the Qa component does not change 
too much, indicating that the variation in the nematic 
ordering mostly occurs in the xz plane, perpendicular to 
the symmetry axis of the particle. Alternatively, Qi, Q2i 
and Qa profiles can be interpreted also in terms of order 
parameters S and P (see Fig. ||, upper panel). When 
the defect line is approached, the uniaxial order param- 
eter S decreases from its temperature-defined bulk value 
and drops even below zero in the defect center. Note 
that there the nematic director, i.e., the eigenvector cor- 
responding to the negative eigenvalue, is directed along 
the long axis of the solid particle. On the other hand, 
the biaxiality - close to zero far enough from the defect 
- increases when the defect line is approached, reaches 
a maximum and, finally, in the very center of the de- 
fect, again drops to a value which is close to zero. The 
characteristic length scales for these variations are of the 
order of a few (w 5) lattice spacings a and agree with 
the estimates for the correlation length in the nematic 
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Qualitatively, molecular ordering close to a dischnation 
line can be summarized as follows. In the very center of 
the defect molecular ordering is uniaxial with S < and 
P — > 0. Far enough from the defect line the nematic liq- 
uid crystal is uniaxial again, however, with S > Q and 
P = 0, as expected in a homogeneous or in a weakly dis- 
torted bulk sample. In the intermediate ring-like region, 
nematic ordering is biaxial with P ^ 0. These conclu- 
sions agree also with results from alignment tensor-based 
phenomenological analyses of topological defects both of 
half-integer 134] and integer strength [^. 



2. MD results 




FIG. 4: MD simulation results: director streamlines of the 
xz cross section of the director field. Rod diameter ar = 
20(70, rod length L = 00. The shading represents the value of 
the density. The director far from the particle is constrained 
along the z axis. A pair of — i line defects forms parallel to 
the particle axis, perpendicular to the director far from the 
particle. 




As has already been mentioned, the configuration with 
two — i dischnation lines is more energetically favorable 
than a single —1 strength dischnation. We noticed this 
performing the molecular dynamics simulation: for all 
studied diameters of the rod {ar/ao = 5 — 20) the — 1 
strength dischnation appears immediately after expand- 
ing the colloid particle in the nematic state. However, 
during the equilibration, it splits into two — ^ dischnation 
lines which then move towards the equatorial plane. The 
evolution dynamics is quite slow, one needs about 10^ 
steps for the — i disclinations to reach the equator. We 
were not able to observe 'escaped radial' configuration, 
probably due to the small size of the colloidal particle. 

A typical director map together with the density map 
is shown in Fig. ^. The bulk density, Pb^o ~ 0.34, is 
slightly different from the number density, pctq = N/V = 



FIG. 5: MD simulation results: order tensor components Qi 
(triangles), Q2 (circles), and Q3 (squares) across the defect. 
In the upper panel we plot the density profile across the defect 
(circles) and avoiding the defect (open squares). The density 
modulation near the particle affects the order parameter vari- 
ation in the core region. 



0.33, because of the volume taken by the rod. In the di- 
rection of the disclinations the density modulation, typ- 
ical for a nematic-wall interface, vanishes due to partial 
melting of the liquid crystal in the dischnation core re- 
gion. This melting damps the influence of the droplet 
surface on the interface region. 

Two — i disclinations are located very close to the 
droplet surface and the director distortion vanishes very 
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quickly in the liquid crystal bulk. The core region ex- 
tends over a few molecular lengths. In MD simulations, 
the pair of defect lines forms perpendicular to the direc- 
tor. Since we use the director constraint algorithm 
the director far from the particle is aligned along the z 
axis, contrary to the situation with MC simulation re- 
sults, where the director is along the box diagonal. The 
director constraint damps the effects of the defect repul- 
sion and effective fluctuation of the director. 

To emphasize the complex structure of the defect core 
we plot the order tensor components Qi, Q2, Q3, after 
diagonalizing the local order tensor, eqn (^, in Fig. 
Qualitatively, the order tensor has the same look as in 
MC simulations: the nematic phase is uniaxial far from 
the core and biaxial in the core region, with variation 
of the biaxiality across the core. However, MD results 
predict a more complicated structure of the core region, 
attributed to the liquid crystal density oscillation near 
the particle surface (see Fig. ^ and Fig. ||, upper panel). 
The variation of the order tensor components is given by 
the superposition of the nematic order variation due to 
the density modulation and intrinsic variation due to the 
presence of the defect. 



B. Torque on the particle 

To measure the torque on the rod of finite length, we 
performed MD simulations in a box with periodic bound- 
ary conditions, applying a global constraint for the direc- 
tor along the z axis ||36| . An independent measurement 
was performed in slab geometry. In the slab geometry, 
the director orientation far from the rod was fixed by the 
confining walls. The walls provided strong homeotropic 
(along the z axis) anchoring of the director. The rod 
was fixed either in the center of the simulation box or 
at some distance z from the bottom wall. > (< 0) 
corresponds to a torque which tends to align the rod per- 
pendicular (parallel) to the director far from the particle, 
no = e^. 

The torque on the rod, calculated using eqn (^, is 
presented in Fig. ^. Results presented in Fig. || indicate 
that the dependence of the torque on the rod tilt angle 
is far from the sin 20, proposed in p^]. Moreover, the 
torque is not equal to zero for Q — , i.e. there is some 
symmetry breaking and the orientation of the rod along 
the director is not even metastable. 

For better understanding, a slice in the yz plane is 
shown in Fig. ^ for different tilt angles of the rod. Fig. |^ 
shows that the director distribution around the rod is not 
axially symmetric. This is the reason for the non-zero 
torque for = 0''. Strong director variations near the 
rod are responsible for the large value of the torque. As 
the rod rotates, the director field becomes less and less 
frustrated. Fig. and finally we have a stable orienta- 
tion of the rod perpendicular to the director. Fig. Qc. 

In principle, the configuration with axial symmetry is 
also possible, when the rod is along the z axis. How- 
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FIG. 6: Torque on the rod vs. rod tilt angle. > 

(< 0) corresponds to a torque which tends to align the rod 
perpendicular (parallel) to the director far from the particle, 
no = fiz. Rod diameter Or = 5(to, rod length L = lOao. 
Squares - slab geometry, with the rod in the middle of the cell. 
Circles - slab geometry, with the center of the particle located 
at the distance z = IScro from the bottom wall. Triangles - 
the depletion force on the particle, when it is near the wall. 
As a guide, the dashed lines correspond to a polynomial fit. 
See also Fig. ^ for explanations. 





FIG. 7: Director streamlines and order parameter maps for 
different tilt angles: a) 61 = 0°; b) f = 45°; c) 6* = 90". A 
side view along the x axis is shown (the rod is tilted in the 
yz plane). Rod diameter a,. = 5(To, rod length L — IOctq. 



ever, we were not able to observe it in our simulations 
even when disordered, isotropic configurations contain- 
ing the colloidal particle were compressed to the ordered, 
nematic state. This method, in principle, gives the lowest 
free-energy configurations in an unbiased way. 

In Fig. ^ we also present the torque and the force on 
the particle near the wall. Besides the torque on the 
particle because of the average molecular orientation of 
the liquid crystal host, there is an effective interaction 
between the particle and the wall. The presence of this 
depletion-like interaction can be understood as a result 
of interactions between the particle and the liquid crys- 
tal molecules which themselves interact with the wall. In 
other words, if the volume close to both the wall and 
the particle overlap, then the host liquid gains accessible 
volume and can increase its entropy. To show that the 
depletion force is responsible for this 'entropic' contribu- 
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tion to the torque, we plot the depletion force (3aoFz{6) 
in Fig. ^. The correlation between the change in the de- 
pletion force and the change in the torque is evident. As 
the rod moves closer to the wall the value of the 'entropic' 
torque increases and can even affect the equilibrium po- 
sition and orientation of the particle. 

For the parameters used in our simulation, the con- 
tribution of the 'entropic' torque only modifies the de- 
pendence of the total torque on the particle tilt angle. 
However, this contribution can dominate, for example, 
when the system is close to the nematic-isotropic tran- 
sition. Then the colloidal particle might have a tilted 
orientation when approaching the wall. Further work on 
these aspects is in progress. 

IV. CONCLUSIONS 

We used molecular dynamics and Monte Carlo tech- 
niques to study a small elongated colloidal particle sus- 
pended in a nematic liquid crystal. Homeotropic bound- 
ary conditions and strong anchoring create a hedgehog 
defect on the particle surface. We have studied the defect 
structure around the particle that cancels this hedgehog 
defect. 

Our simulation results show that, in the case of a 
very long particle, with transverse size much less than 
its length, a configuration with two — i defects is stable. 
An initial configuration with one — 1 strength disclination 
evolves spontaneously into two — disclination lines. For 
a particle with transverse size of the order of its length, 
the — i ring defect encircling the particle was found to be 
stable. The orientation of the ring changes as the particle 
tilts with respect to the director. Using order tensor and 
density maps we are able to resolve the structure of the 
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